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The objective of this grant was to study the transition mechanisms on a flat-plate boundary 
layer interacting with the wake of a bluff body. This is a simplified configuration, sketched in 
Fig. 1, designed to exemplify the phenomena that occur in multi-element airfoils, in which the 
wake of an upstream element impinges on a downstream one. Some experimental data is 
available for this configuration [1 — 3] at various Reynolds numbers. 


Figure 1 Sketch of the geometric configuration. 


The first task carried out was the implementation and validation of the immersed-boundary 
method [4, 5]. This was achieved by performing calculations of the flow over a cylinder at low 
and moderate Reynolds numbers. The low-Reynolds number results are discussed in Ref. 6, 


which is enclosed as Appendix A. The high-Reynolds number results are presented in a paper in 
preparation for the Journal of Fluid Mexhanics. 

We performed calculations of the wake/boundary-layer interaction at two Reynolds numbers, 
Re = 385 and 1155. The first case is discussed in Ref. 6, and a comparison of the two calculations 
is reported in Ref. 7. The simulations indicate that at the lower Reynolds number the boundary 
layer is buffeted by the unsteady Karman vortex street shed by the cylinder. This is shown in 
Fig. 2: long streaky structures appear in the boundary layer in correspondence of the three- 
dimensionalities in the rollers. The fluctuations, however, cannot be self-sustained due to the 
low Reynolds-number, and the flow does not reach a turbulent state within the computational 
domain. In contrast, in the higher Reynolds-number case, boundary-layer fluctuations persist 
after the wake has decayed (due, in part, to the higher values of the local Reynolds number Re 
achieved in this case); some evidence could be observed that a self-sustaining turbulence 
generation cycle was beginning to be established. 



Figure 2 Instantaneous flow’ field for the Re=385 calculation. Iso-surfaces of pressure and streamwise velocity 
fluctuations. 


A third simulation was subsequently carried out at a higher Reynolds number, Re= 3900. This 
calculation gave results similar to those of the /?£■= 1155 case. Turbulence was established at 
fairly low Reynolds number, as a consequence of the high level of the free-stream perturbation. 
Figure 3 shos an instantaneous flow visualization for that case. 

A detailed examination of flow statistics in the transitional and turbulent regions, including 
the evolution of the turbulent kinetic energy (TKE) budget and frequency spectra showed the 
formation and evolution of “turbulent spots”' characteristic of the bypass transition mechanism. 
It was also observed that the turbulent eddies achieved an equilibrium, fully developed turbulent 
states first, as evidenced by the early agreement achieved by the terms in the TKE budget with 
those observed in turbulent flows (Figs. 4 and 5). Once a turbulent Reynolds stress profile had 
been established, the velocity profile began to resemble a turbulent one, first in the inner region 
and later in the outer region of the wall layer. An extensive comparison of the three cases, 
including budgets, mean velocity and Reynolds stress profiles and flow visualization, is in Ref. 
8. 



Figure 3 Instantaneous flow field for the Re=3900 calculation. Iso-surfaces of pressure and streamwise velocity 
fluctuations. 


Copies of Refs. 6 and 7, which discuss the results obtrained, are enclosed. 
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Abstract. The immersed-boundary method can be used to simulate flows 
around complex geometries within a Cartesian grid. This method has been 
used quite extensively in low Reynolds-numb er flows, and is now being 
applied to turbulent flows more frequently. The technique will be discussed, 
and three applications of the method will be presented, with increasing 
complexity, to illustrate the potential and limitations of the method, and 
some of the directions for future work. 


1. Introduction 

The increase in computer speed achieved over the last few years has made 
computational fluid dynamics increasingly useful and widespread as a tool 
to analyze and design flow configuration. Complex geometries, however, 
present an obstacle even to present-day computers, since the use of body- 
fitted meshes (structured or unstructured) significantly increases the cost 
of a calculation in terms of both computational speed and memory require- 
ments. 

An alternative method that may be cost-efficient in many situations is 
the “immersed-boundary” method. This technique is based on the intro- 
duction of body forces distributed throughout the flow that mimic the effect 
that a solid body would have on the fluid. This approach allows the use 
of codes ‘in Cartesian coordinates, which present significant advantages, in 
terms of speed, accuracy and flexibility, over codes that employ body fitted 
grids. 

This idea has been widely used in hsemo-dynamics and bio-fluids en- 
gineering: two- and three-dimensional calculations of the flow in the heart 
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were reported by Peskin [13, 14] and McQueen and Peskin [8, 9]. In these 
calculations the motion of the boundary was determined by the fluid it- 
self, so that the boundary had to be modeled as a set of elements linked 
by springs. In cases in which the boundary motion is known a priori , the 
problem can be significantly simplified. 

Goldstein et al. [4] proposed a feedback forcing mechanism that forces 
the fluid velocity Uj to approach the velocity of the solid boundary. Vy , on 
the boundary itself. Consider the incompressible Navier-Stokes equations: 


duj 

dxj 


= 0, 


dui 


4- 


dt dx 


d \ 1 dp _o . 

(UjUi) — — —— 4 * 4 - 


pdxi 


(1) 

( 2 ) 


Goldstein et al. [4] assigned a force field 

&f J [^1 (r 3j j j t) j £)] dt d~ (3f [v>{ {x$^i , t) ^(^s,:^)]? 

( 3 ) 

where ay and (3f are two negative constants, and x s j are the coordinates 
of the solid surface. The net effect of this force is to tend to annihilate 
the velocity difference U{ — V{. The flow, in fact, responds to the forcing 
as a damped oscillator (see [4] and [3] for an in-depth discussion of this 
issue); the frequency of the oscillator is a lay] 1 / 2 , whereas its damping 
coefficient is cc 0f/\a\ 1 / 2 . This implies that, in order to enforce the no-slip 
condition effectively, ay and (3f must have large magnitudes (larger than 
the highest frequency in the flow), which make the equations stiff. If the 
forcing is advanced in time explicitly, in fact, the maximum allowable CFL 
is o(10" 3 ); even implicit methods only allow the use of CFL numbers of 0.1 
or less. This constraint makes this method impractical for the calculation 
of unsteady (in particular, turbulent) flows. 

Recently, Mohd-Yusof [11] proposed the “direct forcing method,'’ which 
assigns a force field given by 


, d 1 dp o 

— (u jUt ) + -— Ui + 


p dx{ 


Vi ~ m 
!it 


( 4 ) 


(where the dependence on x S) { and t has been omitted). This forcing im- 
poses directly the desired velocity on the immersed boundary, and has the 
advantage (over the feedback forcing method) that it does not require sig- 
nificant reductions in the allowable time-step. It was extensively tested in a 
staggered finite-difference code by Fadlun et al [3], who derived an interpo- 
lation scheme to be used when the boundary does not coincide with a grid 
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Computed velocity ! 
outside the body 

Interpolated velocity 
..at the nearest outside point 1 . 



Figure 1. Interpolation method used to apply the forcing. 

line. Verzicco et al. [15] applied this method to the large-eddy simulation 
of high Reynolds number turbulent flows by calculating the flow inside an 
IC engine. 

In the present paper additional applications of this method will be pre- 
sented, and the potential and limitations of the technique, as well as issues 
that require further study, will be discussed. After a brief review of the 
governing equations and of the numerical scheme used, three test cases will 
be shown: a low Reynolds-number flow over a cylinder in the presence of a 
moving surface (Wannier [16] flow), the flow over a circular cylinder at low 
Reynolds number, and the bypass transition on a flat plate caused by the 
interaction between the boundary layer and the wake of a circular cylinder. 

2. Problem formulation 

Governing the flow are the incompressible Navier-Stokes equations (12). 
The flow solver is a standard 2nd-order accurate method on a staggered 
mesh [1]. The fractional time-step method [2, 6] is used and a 2nd-order 
accurate Adams-Bashforth method is employed for the time advancement. 
A non-reflecting boundary condition [12] is used at the outflow, and periodic 
boundary conditions in the spanwise direction. The inflow and free-stream 
conditions depended on the case studied. 

The direct forcing (4) was used. Since the immersed body does not 
follow a grid line the interpolation method proposed by Fadlun et al. [3] is 
used. The forcing is imposed not at the surface itself, but at the first point 
outside it (see Fig. 1) and the solid body velocity in (4) is replaced with 
the velocity Vi obtained by a linear interpolation between the computed 
fluid velocity outside the body, u t -, and the desired body velocity Vi. This 
method has several advantages: first, it has been shown to be fully second- 
order accurate in time [3]; therefore, it is consistent with the second-order 
differencing scheme used by the solver; secondly, since it assumes that the 
velocity profile is linear near the body, it implies homogeneous Neumann 
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Figure 2. Wannier flow test case (Wannier 1950). (a) Computational domain and com- 
puted streamlines; (b) L\ and L 2 norms of the error, e, for u and v velocity components. 
N is the total number of grid points. 


boundary conditions for the pressure (see the Appendix in [3] for a full 
discussion of this issue). This last feature is very important in the framework 
of fractional time-step methods, since it implies that the corrector step does 
not result in a modification of the body velocity imposed through the forcing 
in the Helmholtz step. On the other hand, the assumption that the velocity 
profile is linear over the first two layers of cells outside the immersed body 
requires the use of a very fine mesh in the vicinity of the body. 

3. Results and discussion 

3.1. WANNIER FLOW 

A straightforward way to verify the accuracy of the proposed methodology 
is to compute a flow containing a curved immersed boundary for which 
an analytical solution exists. The case considered here is the Stokes flow 
around a cylinder in the vicinity of a moving wall (see Fig. 2). An analytical 
solution for this case was derived by Wannier [16]. The streamlines for this 
Sow axe shown in Fig. 2a. Three computations on gradually finer uniform 
grids (32x32, 64x64, and 128 x 128) were conducted. The L\ and Lo norms 
of the error (the difference between the computed and analytical solution) 
axe shown in Fig. 2b as a function of the total number of points N. The 
error decreases with a —2 slope indicating that the proposed methodology 
is second order accurate. 

3.2. FLOW OVER A CIRCULAR CYLINDER 

The next test case examined is the flow over a circular cylinder at Reo = 
U^Djv = 300 (where D is the cylinder diameter and Uqq the tree-stream 
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Figure 3. Flow over a circular cylinder, R&d = 300. Velocity statistics, (a) U, (b) W, 
(c) (uV), (d) {n't o'). Reference data from Ref. [10]. 


velocity). Two calculations will be compared: a 2D one that used 400x200 
points in the arz—plane, and a 3D one, with the same mesh in the zz— plane, 
and 48 points in y. The computational domain was 60 x 2 tt x 30, and the 
cylinder center was at x c — 10, z c — 15 (all lengths are made dimensionless 
by D, all velocities by Uqq). The grid was stretched both in the x— and 
z— directions; the last 1/6 of the domain (which required only 10 grid points 
in x) formed a sponge region, used to minimize the upstream propagation of 
disturbances due to the convective outflow conditions. A uniform velocity 
profile was imposed at the inlet, and slip- wall conditions were applied at 
z = 0 and z — 30. 


The velocity statistics are shown in Fig. 3. Here and in the following the 
angle brackets denote averaging in time and in the spanwise direction. The 
3D calculation is in very good agreement with the reference data by Mittal 
and Balachandar [10]. At this Reynolds number, three-dimensionality is 
observed in the wake, which is evidenced in the visualization in Fig. 4, 
which shows iso-surfaces of the second invariant of the velocity- gradient 
tensor, 


1 dui duj 

2 dxj dxi 


-(SijSij-sUjnij), 


( 5 ) 
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(where is the anti-symmetric part of the velocity gradient tensor). The 
condition Q > 0 identifies effectively the regions of coherent vorticitv 5 . 

Figure 4 shows the formation of an instability on the initially 2D rollers, 
and the presence of quasi-streamwise rib vortices joining the rollers. The 
magnitude of the spanwise Reynolds stresses (uV). in this calculation, how- 
ever. remained significantly smaller than the other two normal components, 
which allowed the 2D calculation to give reasonable results. 

The effects of the grid resolution near the obstacle are shown in Fig. 5. 
The cell Reynolds number is defined as 


Re c 


(Ax 2 -f Az 2 ) ( u 2 

V 



(S) 


where u and w are the instantaneous velocities, and Ax and Az the grid 
spacings. If the mesh is insufficiently fine {Re c > 30), some oscillations 
can be observed that initiate along lines at ±45° on the cylinder (they 
are especially visible in the w contours. Fig. 5b). Refining the grid, thus 
reducing Re c reduces the size of this oscillation (Figs. 7a and b). Two- 
dimensional interpolation schemes that use both the points indicated by 
the diamonds and those indicated by squares in Fig. 6 to determine V x 
have also been found (Verzicco, private communication, 2001) to reduce 
the amplitude of these oscillations. 
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Figure 5. Flow over a circular cylinder, ReD = 300. Coarse (200 x 100) 2D calculation, 
(a) Contours of the cell Reynolds number; (b) contours of w. 



Figure 6. One-dimensional vs. two-dimensional interpolation stencils. 


3.3. WAKE/BOUNDARY-LAYER INTERACTION 

Wakes interact with laminar boundary layers in many applications of en- 
gineering interest, for example on the leading edge of multi-component 
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Figure 7. Flow over a circular cylinder. Reo = 300. rine (400 x 200) 2D calculation, 
(a) Contours of the cell Reynolds number; (b) contours of w. 



Figure 3. Sketch of che wake/boundary- layer configuration. 


airfoils, or inside turbo-machinery. The interaction of the turbulent eddies 
present in the wake with the boundary layer itself r hen becomes a primary 
driver of the transition process in the boundary layer itself, and may lead 
to transition to turbulence at fairly low Reynolds numbers. 

The configuration examined in the present study is sketched in Fig. 3. 
A circular cylinder, with its axis normal to the stream is placed above a 
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flat plate. The cylinder center is at x c = 10, z c = 3.2, immediately above 
the leading-edge of the plate, which was also at x c = 10. As in the previous 
case, distances are normalized by D , velocities by U^. The computational 
domain was 60 x 27r x 20. As for the cylinder calculation, the grid was 
stretched both in the x— and z- direct ions and a sponge region was used. 
The Reynolds number based on cylinder diameter was 385. The configu- 
ration corresponds to Case 1 in the experimental paper by Kyriakides et 
al. [7], who observed significant velocity fluctuations in the boundary layer, 
starting from a point approximately six diameters downstream of the cylin- 
der. These fluctuations are generated by the large-scale convective motion 
of the vortices, and do not die down after the wake has weakened, but de- 
velop into a turbulent boundary layer despite the fact that the Reynolds 
number is very low. 

The distribution of the streamwise Reynolds stress, (uV), as a function 
of x is shown in Fig. 9. A sudden increase of (uV) can be observed to begin 
at x — 8, indicating the beginning of transition. This result is consistent 
with the observations of Kyriakides et al [7], who defined the onset of 
transition as “the x ~ location where the velocity signal at the same height 
above the plate loses its sinusoidal character”, and found that transition 
occurs at x = 7.4. 

The velocity profiles, shown in Fig. 10a at several locations, initially 
resemble a Blasius profile merging into a wake near the cylinder. As the 
wake widens and interacts with the boundary layer, a logarithmic layer 
begins to establish itself, indicative of transition towards turbulent flow. 
This transitional behavior is also observed in the trace of the Reynolds 
stress tensor, q 2 (equal to twice the turbulent kinetic energy), which in 
the latter sections establishes a turbulent-like distribution, with a peak of 
magnitude 7 — 8 r w at z + ~ 10 — 12. It should be noted that this quasi- 
turbulent state is achieved at very low Reynolds number: the boundary- 



Figure 9. Wake/ boundary- layer interaction, (uv!) distribution at z = 0.13. 
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Figure 10. Wake /boundary- layer interaction. Turbulent statistics. Top: mean velocity 
profile; bottom: q 2 = 

layer thickness 5 (defined as the distance above the plate at which the first 
maximum of the velocity profile occurs) is approximately 50-70 wall units. 

A visualization of the flow is shown in Fig. 11. The structure of the 
cylinder wake is similar to that highlighted in Fig. 4, with strong spanwise 
rollers that exhibit 3D instabilities and eventually break up, and smaller 
quasi-streamwise vortices in the braid region. The contours of streamwise 
velocity fluctuation a' on a plane parallel to the wall show significant; levels 
of fluctuations, especially near kinks in the rollers belonging to the lower 
row. Quasi-streamwise streaks are formed around x = 15, whose spacing is 
approximately 100 wall units. 

4. Conclusions 

The immersed-boundary technique has been presented and discussed. Il- 
lustrative results from three simulations indicate the potential of this tech- 
nique, which allows the calculation of flows around complex geometries 
without requiring a body-fitted grid. 
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Figure 11. Wake /boundary-layer interaction. Isosurfaces of Q = 0.4 and contours of u 
in the 3 = 0.13 plane. Top: perspective view; bottom: view from above. 


If appropriate interpolation methods are used when the body does not 
coincide with a grid line, the method is second-order accurate. However, 
some care must be taken in the discretization of the Sow held, especially in 
the vicinity of the body. We have observed the development of numerical os- 
cillations when the cell Reynolds-number exceeded values of approximately 
30. These oscillations did not grow or change position in time, and their 
effect on the Sow field downstream of the obstacle was limited in the cases 
studied. This w^as observed, for instance, in calculations of the Sow around 
the cylinder at Rejg = 3500 were carried out in which f he grid could not be 
sufficiently renned to satisfy the cell- Reynolds-number requirement. How- 
ever, it is not known whether these oscillations might give rise to instabili- 
ties in other geometries, or at higher Reynolds numbers. The development 
of more accurate, multi-dimensional interpolation schemes might be bene- 
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ficial in this respect. The use of multi-block methods, or embedded grids, 
could also alleviate this problem. 

If these numerical schemes can be overcome, the immersed ccundary 
method appears to be a useful tool for the simulation of flows in complex 
geometries at moderate or high Reynolds numbers. This is confirmed by 
the increasing number of studies using this method that are appearing in 
the literature. 
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Direct and large-eddy simulations of the interaction between the wake of a circular 
cylinder and a flat-plate boundary iayer are conducted. Two Reynolds numbers are ex- 
amined. The simulations indicate that at the lower Reynolds number the boundary layer 
is buffeted by the unsteady Karmdn vortex street shed by the cylinder. The fluctuations, 
however, cannot be self-sustained due to the low Reynolds-number, and the flow does 
not reach a turbulent state within the computational domain. In contrast, in the higher 
Reynolds-number case, boundary-layer fluctuations persist after the wake has decayed 
(due, in part, to the higher values of the local Reynolds number Re$ achieved in this 
case); some evidence could be observed that a self-sustaining turbulence generation cycle 
was beginning to be established. 


Introduction 

High-lift systems have a significant impact on the 
overall cost and safety of aircraft. According to Mered- 
ith, 1 a 1% improvement in the maximum lift coefficient 
(or lift-to-drag ratio) could translate into an increased 
payload of 14 to 22 passengers on a large twin-engine 
transport. An optimal aerodynamic design of a multi- 
air foil high- lift configuration requires careful consider- 
ation of both inviscid and viscous flow phenomena. In 
particular, laminar-to-turbulent transition is a crucial 
issue for ground-to-flight scaling of high-lift flow fields. 

Either of the familiar transition mechanisms for a 
single-element configuration are also relevant to tran- 
sition over a multi-airfoil configuration: those asso- 
ciated with streamwise instabilities in the form of 
Tollmien-Schlichting or Rayleigh modes, cross-flow 
and attachment line instabilities, and leading-edge 
contamination — a form of subcritical (he., bypass) 
transition. However, a unique transition mecha- 
nism in the case of multi-element flow-fields involves 
the boundary-layer contamination due to unsteady 
wake(s) from the preceding element(s) and/or addi- 
tional forms of vortical disturbances originating from 
the separated, cove flow underneath an upstream el- 
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ement (Fig. 1). While the “single-element” class of 
transition mechanisms have been widely studied in the 
literature, the wake-contamination issue has received 
little scrutiny thus far and is the focus of this paper. 

Interactions between turbulent wakes and boundary 
layers have been the subject of much study. Most of 
the investigations, however, concentrated on the un- 
steady case that occurs, for instance, in turbomachin- 
ery, when the wake of an upstream blade impinges on 
a downstream blade; the unsteadiness of the impinge- 
ment region plays an important role in the dynamics 
of the flow. Fewer studies can be found of the steady 
case. Squire 2 summarizes much of the investigations 
conducted prior to 1989. Particularly important is the 
work by Zhou and Squire 3 who examined the interac- 
tion between the wake of am airfoil and a flat plate. 
They found that, a region in which the wake and 
boundary layer are separated by a potential core is fol- 
lowed by a merging zone, in which the velocity profile 
in the outer layer is substantially different from that 
in a regular flat-plate boundary layer. They also ob- 
served that the position of zero Reynolds shear stress 



Fig. 1 Sketch of the wake/boundary-layer inter- 
actions on an airfoil. 
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and that of zero mean-velocity gradient do not coin- 
cide, an important issue for eddy-viscosity turbulent 
models. 

Recently, Kyriakides and co-workers 4 ' 6 performed 
a series of experiments involving the interaction be- 
tween the wake of a large- aspect-ratio circular cylinder 
and a flat-plate boundary layer. These experiments 
spanned a range of Reynolds numbers and cylinder 
diameters, including cases in which the wake still con- 
tained coherent eddies in the interaction region (i.e., 
the so called case of “strong interaction”) and oth- 
ers in which the interaction took place sufficiently far 
downstream of the cylinder, such that the spanwise 
rollers had decayed significantly (i.e., “weak interac- 
tion” cases) . While they provided measurements for a 
range of Reynolds numbers, Kyriakides et al . A ~ 6 did 
not clarify how transition could occur (and turbulence 
could be self- sustained) at the lowest Reynolds num- 
bers examined. 

Due to the tremendous challenges inherent in a high- 
fidelity simulation of the wake-contamination problem, 
we also choose to study a simple, building-block prob- 
lem, namely that of interaction between a flat- plate 
boundary layer and the wake of a circular cylinder 
that is placed above the flat plate. This model prob- 
lem provides a reasonable balance between simplicity 
(two-dimensional geometry involving a combination of 
two canonical flow fields) and the numerical challenges 
of a complex-geometry simulation including the wake 
generator and a boundary layer. Simultaneously, it al- 
lows one to investigate a range of issues, such as the 
effect of Reynolds number and the role of the spac- 
ing between the cylinder and plate (which allows one 
to traverse the continuous range from strong wake 
interactions to weak interactions that are analogous 
to conventional free-stream turbulence). Effects of 
flow three-dimensionality and/or transverse pressure 
gradients (which lead to curved wakes) may also be 
captured, if necessary, by modifying the free-stream 
boundary conditions. Most significantly, a limited 
amount of experimental data is available for a range 
of flow and geometry parameters . 4 " 6 

In the present paper we report the initial re- 
sults from a series of numerical experiments de- 
signed to study the physical phenomena underlying 
wake/boundary-layer interactions of this type. The 
first stage of this investigation involves the study of 
the strong interaction case, with the specific aims to 
clarify some issues left unresolved by the experimen- 
tal studies, and to document the flow field in detail. 
In particular, we will focus on the issue of transi- 
tion mechanisms at low Reynolds number, which was 
only partially addressed in the experiment. Thus, we 
choose to examine low Reynolds-number cases first, 
before attempting the more challenging simulations at 
higher Reynolds numbers. Future work will focus on 
the weak interaction case, as well as other shapes of 


wake generators. 

Herein, we analyze direct and large-eddy simulations 
of two cases studied by Kyriakides et a /., 4 having the 
same geometry but different Reynolds numbers, which 
lead to different responses of the boundary layer to the 
unsteady wake in the free stream. Despite the simplic- 
ity of the geometrical configuration, conducting DNS 
and LES of the above problem is a challenge for most 
numerical methods available today. In the present 
study we model the effect of the cylinder using an “im- 
mersed boundary” formulation . 11 This approach en- 
ables one to use codes in Cartesian coordinates, which 
present significant advantages, in terms of speed, accu- 
racy and flexibility, over codes that employ body-fitted 
grids. 

The outline of the paper is as follows: the numeri- 
cal formulation of the problem including the geometry 
of the cylinder-plate configuration are described first. 
Then, the numerical methodology will be presented. 
We will then present validation of the numerical tool 
in the context of two building-block problems rele- 
vant to the flow configuration of interest; results of the 
simulations will be compared with numerical and ex- 
perimental data. Results for the wake/ boundary-layer 
interaction are discussed next. Finally, conclusions 
and recommendations for future work are presented. 


Problem formulation 

In this work we present both DNS and LES re- 
sults. In the DNS case, the Navier-Stokes equations 
are solved: 


duj 

dxj 


- 0 , 


a) 



9 ( \ - 1 dP 

dxj UjUl ' pdxi 


+ i/V 2 Ui + fi. (2) 


For the LES, we use the filtered equations of conser- 
vation of mass and momentum 


dxj ’ 

dUi d _ _ 

+ (u jUi ) - 

dru , -7 


dt 


dp 

dx x 


dxj 





(3) 

(4) 


(w’here the over-bar denotes filtered variables and the 
effect of the subgrid scales appears through the SGS 
stresses r X j = ujuj — UiUj)] fi and f i are the body 
forces used in the immersed-boundary method to en- 
force the no-slip conditions on the body. 

For both DNS and LES, the equations of motion 
are solved numerically using a second-order accurate 
finite-difference method on a staggered grid. In all 
cases the grid is uniform in the spanwise direction 
y 7 stretched in the streamwise and normal ones (x 
and z, respectively) to allow accurate resolution of 
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Fig. 2 Sketch of the computational configuration 
and boundary conditions. The drawing is not to 
scale. 

the boundary layer and of the shear layers emanat- 
ing from the cylinder. The discretized equations are 
integrated in time using an explicit fractional time- 
step method, 7,8 in which ail terms axe advanced in 
time using the Adams-Bashforth method. The Pois- 
son equation is then solved using a direct solver, and 
the velocity is corrected to make the field solenoidal. 
For the LES, the SGS stresses are parameterized using 
the Lagrangian Dynamic eddy- viscosity model, 9 which 
has been shown to give accurate results in transitional 
flows. 

The effect of the cylinder, which does not coincide 
with the grid lines, is introduced using an “immersed 
boundary” formulation. In our implementation the 
enforcement of non-slip boundary conditions on the 
cylinder surface is performed using a method similar 
to the the “discrete forcing” approach introduced by 
Mohd-Yusof 10 and Fadlun et al . 11 With this approach, 
the solution is reconstructed at the grid nodes near the 
boundary to ensure satisfaction of the no-slip condi- 
tions. While Mohd-Yusof 10 and Fadlun et al n used a 
simple one-dimensional reconstruction scheme, in the 
present study we employ a multidimensional proce- 
dure proposed by Balaras, 13 in which the solution is 
reconstructed along the line normal to the interface. 

The code is parallelized using MPI message-passing 
routines; the parallelization is performed by domain- 
decomposition process. Except for the Poisson solver, 
the computational domain is divided into a number 
subdomains in the z directions, each of which is as- 
signed to one processor. In the Poisson solver, after 
the velocity divergence is Fourier- transformed in the 
span-vise directions, each processor is responsible for a 
set of spanwise modes. 

We applied the following boundary conditions (see 

Fig- 2): 

1. At the inlet, a uniform velocity in the x direction 
was assigned. 

2. At the outlet, radiative boundary conditions 17 
were given. Furthermore, over the last 30% of 
the domain (which contains only 5% of the grid 


points) the mesh was stretched very significantly, 
to generate a buffer layer in which the equations 
are effectively parabolized because of the large as- 
pect ratio of the cells. 

3. In the spanwise direction, y, periodic conditions 
were used. 

4. In the ffeestream, we used a slip-wall boundary 
condition. 

5. On the cylinder surface, no-slip conditions were 
enforced using the methodology described in 
Balaras. 13 The cylinder was placed at a height of 
3.5D above the nominal leading-edge of the plate. 

6. On the z — 0 surface, we used free-slip condi- 
tions (w — 0, du/dz — dv/dz — 0) ahead of 
the plate, and no-slip conditions (u< = 0) on the 
plate. A sharp transition between free-slip and 
no-slip was found to generate numerical oscilla- 
tions; therefore, a hyperbolic tangent profile was 
used to merge free-slip and no-slip conditions in 
a smooth fashion. As a result, the origin of the 
boundary layer on the plate cannot be determined 
precisely; in the present calculations, the virtual 
origin of the boundary layer appears to be ap- 
proximately l.oD upstream of the actual leading 
edge. 

Two calculations were carried out, corresponding to 
cases 1 and 2 of Kyriakides et al . 4 The geometry for 
the two cases was the same, but the free-stream ve- 
locity was varied to obtain Reynolds numbers Reo 
equal to 385 and 1,155, respectively. The definition 
of the Reynolds number is based on the free-stream 
velocity and cylinder diameter. The extent of the com- 
putational domain in the streamwise and wall-normal 
directions was 100D and 30D, as shown in Fig. 1. The 
spanwise width of the domain was 2i tD wide, a length 
sufficient to include several rib vortices, and to resolve 
the three-dimensional structures in the cylinder -wake. 

Code validation 

The code used in this paper has been exten- 
sively validated for a variety of turbulent 14,10 and 
re-laminarizing 18 flows. The “immersed boundary” 
formulation which is used to enforce non-slip boundary 
conditions on the cylinders surface, has also been ex- 
tensively tested in the context of the present code. 13 In 
particular, a series of three-dimensional computations 
were conducted of the flow around a circular cylin- 
der at Reo = 300. These involved a total number of 
grid points ranging from 1.5 x 10 6 to 4.7 x 10 6 . The 
agreement of the results with reference computations 
using boundary fitted coordinates demonstrated the 
accuracy of the method. Details on the above compu- 
tations can be found in Balaras. 13 
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Fig. 3 St Teamwise development of the skin- friction 
coefficient Cf. 



Fig. 4 Profiles of the streamwise turbulence in- 
tensity, Urms • 

The grid required for accurate simulations of bypass 
transition needs to be determined carefully. To this 
end, we performed simulations of bypass transition on 
a flat plate due to free-stream turbulence. We used the 
setup adopted by Voke and Yang. 18 The calculation 
was started at the streamwise location Re x ~ 6, 000; 
the inflow condition consisted of a Blasius boundary 
layer, on which turbulent fluctuations were added. 
The fluctuations were obtained from a separate calcu- 
lation of spatially developing homogeneous isotropic 
turbulence. To avoid introducing excessively large 
fluctuations near the wall, the turbulent fluctuations 
were windowed in such a way that they vanished at 
25 0 /3 (where 5 0 is the inflow boundary-layer thick- 
ness). Two cases were run, with a 3% and a 6% 
freestream turbulence amplitude, to match the calcu- 
lations 13 and the experimental data. 19 The descrip- 
tion of the inflow conditions both for the experiment 
and the LES calculations was not, however, detailed 
enough to allow us to reproduce exactly the experi- 
mental setup. 

We performed DNS of both cases. For the 3% case, 
we 850 x 96 x 96 points to simulate a domain that 
had dimensions 714<5 0 x 32 5 0 x 4 S5 0 ; For the 6% case, 
550 points were used in the streamwise direction to 
discretize a box that was 4765 0 long. The mesh was 


uniform in the spanwise direction, and stretched in the 
x and z directions. In the transition region, the grid 
size was approximately Ax = 0.25, Ay = 0.085; at 
least 40 points were used to discretize the local bound- 
ary layer thickness, 5. For the LES of the 6% case, 
300 x 64 x 96 points were used, resulting in coarser 
resolution in x and y: Ax = 0.45, Ay = 0.15. 

The streamwise development of the skin-friction co- 
efficient Cf — 2r w / pUl is shown in Fig. 3. One can 
observe a higher skin-friction coefficient in the lami- 
nar region, due to the perturbation introduced at the 
inflow, which generates a pressure disturbance that 
propagates well inside the boundary layer and dis- 
torts the mean- velocity profile in the near- wall region. 
It was verified that, in the absence of inflow distur- 
bances, the Blasius profile was recovered. The slope 
of the skin-friction coefficient distribution is, however 
correct. The transition begins at the correct location 
for both values of the freestream turbulence intensity, 
and the extent of the transition region is also predicted 
with reasonable accuracy. The LES calculation pre- 
dicts the a slightly delayed beginning of the rise of 
Cf, and a somewhat more rapid transition; the results 
are, however, in reasonable agreement with the DNS 
and the experiments. In Fig. 4 the distributions or the 
r.m.s. turbulence intensities at Re x — 290, 000 (for the 
3% freestream turbulence case) are compared with the 
experimental data of Roach and Brierley. 19 The agree- 
ment is again satisfactory, given the uncertainty in the 
inflow condition generation. 

Finally, a grid refinement study was carried out on 
the real geometry. Two grids were tested, a coarse one 
that used 576 x 48 x 192 grid points and a medium one 
with 864 x 72 x 288 points. Calculations on a finer grid 
with 1056 x 128 x 384 points are currently underway. 
The coarse grid was similar to the one that gave grid- 
converged results in the cylinder calculation. Figure 5 
shows the mean velocity profile and the trace of the 
Reynolds stresses q 2 = (u^u-) (throughout this paper 
(■) represents a long-time average, and /' = / - (/) 
is the fluctuating component of /) at three locations 
downstream of the cylinder. Very little difference can 
be observed between the coarse and medium grids for 
the mean velocity profiles; the coarse resolution of the 
wake and boundary layer results in some difference in 
q 2 between the coarse and medium meshes, but the 
agreement is altogether satisfactory. In the following, 
the medium grid results will be presented. 

Notice that the resolution of the boundary layer in 
the medium grid compares favorably with the bypass 
transition case mentioned above. In the transition re- 
gion, the grid spacing was Ax = 0.075, Ay — 0.065 
(where 5 is the local boundary layer thickness). At 
least 80 points were used to discretize the boundary- 
layer thickness. 
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Fig* 6 Mean velocity contours. Ren = 385; Ren = 1, 155. 
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Fig. 5 (a) Mean velocity profiles and (b) profiles 

of q 2 , Rsd = 385. 576 x 48 x 192 grid; 864 x 

72 x 288; X Blasius solution. 
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x/D 

Fig. 7 Streamwise development of: (a) the ve- 

locity at the edge of the boundary layer; (b) the 
acceleration parameter A'. 

Results and discussion 

The mean- velocity contours from the 

wake/boundary-layer simulations at Ren — 385 

and Reo = 1, 155 are shown in Fig. 6. As expected, 
both the wake and the boundary layer are thinner 
at the higher Reynolds number. No merging is 
observed within the computational domain. Due to 
the presence of the cylinder, the boundary layer near 
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Fig. 8 Streamwise development of statistical quan- 
tities. (a) skin friction coefficient C/; (b) Reynolds 
number based on the momentum thickness, Res', (c) 
peak streamwise and spanwise Reynolds stresses in 
the boundary layer. 

the leading edge first encounters an accelerating outer 
flow, and then a decelerating one. The corresponding 
acceleration parameter K = (i'/U 00 )(dU 00 /dx) (where 
Uco is the velocity at the edge of the boundary layer) 
is initially large and positive. The region of large 
positive values of JC is, however, limited to the region 
x/D < 1 (Fig. 7); over most of the flat plate the 
acceleration parameter is below the value that may 
result in re-laminarization (K ~ 3 x 10~ 6 ). In fact, 
for x/D > 22 the boundary layer is subjected to a 
mild adverse pressure gradient that is expected to 
favor transition to turbulence. 

The streamwise development of several statistical 
quantities is shown in Fig. 8. The skin friction coeffi- 
cient Cf (Fig. 8a) from both calculations is very close 
to the value corresponding to the Blasius boundary 
layer. In the high Reynoids-number case we observe a 
significant increase of Cf (towards the end of the com- 
putational domain); this issue will be discussed further 
below. The Reynolds number based on the momentum 
thickness, Res, is shown in Fig. 8b. For both calcu- 
lations Res remains very low. As mentioned before, 
the Reg values achieved in the Rec = 385 case are 
unquestionably below the minimum values at which 
self-sustained turbulence is typically observed in flat- 



Strouhal 


Fig. 9 Power spectra of the streamwise velocity 
fluctuations near the wall, (a) Data of Kyriakides 
et ai.; 4 (b) present calculation, Reo = 385. 

plate boundary layers. The Reg values from the high 
Reynoids-number calculation towards the end of the 
computational domain are perhaps marginal in this 
regard. 

The power spectra of the streamwise velocity at 
z/D — 0.15 are compared in Fig. 9 with the exper- 
imental data of Kyriakides et al 4 The experimental 
data was taken from Fig. 7 of the paper cited, which 
did not specify the distance from the wall at which 
the data was measured, nor the scale of the ordinate; 
the frequency has been rescaled to yield a dimension- 
less Strouhal number. Our calculations are in general 
agreement with the experimental observations; the 
shedding frequency is 0.22 (in the experiment a value 
of 0.21 was observed) and we observe a filling-up of the 
spectrum at both lower and higher frequencies with 
distance downstream of the cylinder. We do not ob- 
serve, however, the harmonic peak measured in the 
experiment at St = 0.42. It is unclear whether this 
discrepancy is due co experimental or numerical er- 
rors, to insufficient sample convergence in the DNS, or 
to a difference in the location of the measurements. 

Kyriakides et al 4 observed that time-traces of the 
velocity at this height had a sinusoidal behavior up 
to x/D ~ 7.3 and 2.7 for the low and high Reynoids- 
number, respectively; they identified transition with 
the loss of the sinusoidal behavior. The time histories 
of the velocity fluctuations (not shown) show depar- 
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Fig. 10 Mean velocity profiles at various locations. 

Re D = 385; R&d = 1, 155; x Biasius profile. 

Each profile is shifted by 1 unit for clarity. 

ture from the sinusoidal behavior at similar locations. 
We believe, however, that a quantitative definition 
of the onset of transition is more desirable than the 
subjective one adopted by Kyriakides et aJ. 4 One pos- 
sibility is to define the location of the transition point 
in terms of some threshold value of the turbulent fluc- 
tuations. Figure 8c shows the developments of the 
peak (u'u 1 ) and (vV) Reynolds stresses, normalized 
by the local wall stress r w . The streamwise Reynolds 
stresses reach values comparable to those expected in 
turbulent flows in the Rep = 1,155 case, while they 
remain at about half of this level when Reo = 385. 

An important measure of transition is the onset of 
three-dimensionality. We observe that the spanwise 
Reynolds stress (vV) begins growing earlier for Rep = 
385 than for Reo — 1, 155. This is a consequence of 
the faster spreading of the cylinder wake at the low 
Reynolds-number. It will be shown later that the first 
onset of three-dimensionality in the boundary layer is 
caused by the passage of the vortices shed in the wake 
of the cylinder. The onset of three-dimensionality in 
the boundary layer is very well correlated with the 
onset of three-dimensionality in the wake. 

Figures 10 and 11 show the profiles of the mean ve- 
locity and q 2 at several locations downstream of the 
cylinder. We observe a clear separation between the 
boundary layer and the wake near the leading edge, 
which disappears as one moves downstream. Up to 
x/D — 20 the velocity profiles agree very well with 
the Biasius solution at both Reynolds numbers. In the 
low Reynolds-number case we observe a mild growth 
of fluctuations in the boundary layer (evidenced both 
by the streamwise development of (uV), Fig. 8, and 
by the profiles of q 2 , Fig. 11). These fluctuations, how- 
ever, appear to be due entirely to the advcction caused 
by the vortices shed by the cylinder, which, as they are 



Fig. 11 Profiles of q 2 at various locations. 

Reo = 385; Rgd = 1, 155. Each profile is shifted 

by 0.05 units for clarity. 



Fig. 12 Profiles of q 2 in wall units at various lo- 
cations. Ren = 385; Reo - 1, 155. Each 

profile is shifted by 5 units for clarity. 

convected downstream, induce ejections of low-speed 
fluid from the near- wall region. These velocity fluctu- 
ations, in fact, do not result in the levels of Reynolds 
shear stress ( u'w ') typical of the turbulent flat-plate 
boundary layer (see below). 

The high-Reynolds-number calculation shows a dif- 
ferent behavior. Significant differences appear in the 
velocity profile at the downstream locations. At the 
last location shown {x/D = 55), Re x ~ 58,000 and 
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Fig. 13 Profiles of {u w) in wall units at various 

locations. Reo = 385; Reo = 1, 155. Each 

profile is shifted by 1 unit for clarity. 

Res* — 340. Although the flow has not achieved a fully 
developed turbulent state yet, there are several indica- 
tions that transition to turbulence is taking place. The 
turbulent kinetic energy (Fig. 11), for instance, shows 
a near- wall peak that is increasing and moving towards 
the wall. The value of this peak and its location in 
wall units are those expected in a turbulent flat-plate 
boundary layer (Fig. 12). The Reynolds shear stress 
{u , w t ) (Fig. 13) exhibits a similar behavior: its peak 
value and distribution for the high-Reynolds-number 
case are similar to those expected in turbulent flat- 
plate flows, while in the Reo — 385 case the peak is 
significantly lower. 

The Reynolds stresses appear to move towards a 
turbulent state more rapidly than the mean velocity 
profile; the shape factor H = 5*/ 9 remains close to the 
laminar value of 2.6 throughout the computational do- 
main, and the mean velocity profile has not achieved a 
logarithmic profile by the end of the computational do- 
main (Fig. 14). However, significant levels of turbulent 
fluctuations (r.m.s. levels between 1 and 10% of the 
freestream velocity) are observed within the bound- 
ary layer. At the downstream stations, the freestream 
turbulence at the edge of the boundary layer is very 
significant, between 3 and 5% of the velocity at the 
boundary-layer edge (for the low and high-Reynolds 
number cases, respectively), levels that are able to 
trigger bypass transition in boundary layers. Com- 
putations that use a longer domain should be carried 
out to determine the exact route to turbulence in the 
present configuration. 
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Fig. 14 Mean velocity profiles in wall units at 

various locations. Reo = 385; Reo = 1,155. 

Each profile is shifted by 10 units for clarity. 

Figures 15 and 16 show instantaneous contours of 
the streamwise velocity fluctuations, and iso-surfaces 
of the second invariant of the velocity gradient tensor, 

« = = (5) 

where S TJ is the strain- rate and H x; - the rotation-rate 
tensor. In regions where Q > 0, the vorticity is due to 
rotational motions, rather than to shear. 

In the low- Reynolds- number case (Fig. 15) the cylin- 
der wake is very well-defined; staggered row's of quasi- 
2D vortices are shed, which are joined by clearly ob- 
servable rib vortices. Three-dimensionalities quickly 
develop. There is a high correlation between the ap- 
pearance of three-dimensionalities in the rollers, and 
the appearance of low- and high-speed streaks in the 
boundary layer (a very clear example can be observed 
around x/D — 20). This feature is connected to the 
increase of the spanwise Reynolds stress ( v'v ') men- 
tioned above. These streaks are convected with the 
freestream velocity, confirming the observation that 
they are induced by the outer flow: in turbulent flows, 
the typical near- wall streaks, which are generated in 
the inner layer, are convected at speeds that are 60- 
70% of the freestream velocity. 

Another effect of the shed vortices on the near-wall 
region can be observed in Fig. 15c. Corresponding 
to the vortex street, one can observe alternating re- 
gions of high- and low-speed fluctuations, due to the 
advection from the upper and lower rows of vortices, 
respectively. These regions correspond to Q2 and Q4 
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events, 20 and are responsible for regions of significant 
correlation u'w r . As the rollers become weaker, no new 
streaks are generated, and the flow does not appear to 
evolve into a fully turbulent state, the main distur- 
bances being due to the streaks that are convected 
downstream, and weaken progressively. 

At the high Reynolds-number (Fig. 16) a some- 
what different picture emerges. Smaller-scale vortices 
can be observed in the wake, consistent with the in- 
creased Reynolds number. Also, the wake breaks down 
more rapidly. Well-defined tollers and rib vortices can, 
nonetheless, be observed. Some low-speed regions in 
the boundary layer can still be observed, correspond- 
ing to the onset of three-dimensionality in the rollers; 
unlike the low-Reynolds-number case, however, these 
streaks do not weaken significantly as they are con- 
verted downstream. The decreased viscous dissipation 
due to the increased Reynolds number certainly plays a 
role in this, but a self-sustaining turbulence generation 
cycle may also be beginning to appear. This conjecture 
is supported by the appearance of quasi- stream wise 
vortices in the near- wall region for x/D > 50. 

Conclusions 

The work presented in this paper had a two-fold 
objective, namely, to establish the accuracy of the 
numerical code for wake-boundary-layer interaction 
problems and to clarify unanswered questions regard- 
ing prior experiments 4 related to such interactions. We 
believe that both of these goals have been nearly met 
through the work described here. In addition to code 
validation for building-block problems involving the 
isolated problems of cylinder wake and bypass tran- 
sition in a flat-plate boundary layer, we have carried 
out large-eddy and direct simulations to understand 
the details of the strong interaction between a cylinder 
wake and a flat-plate boundary layer at two Reynolds 
numbers (Re d ~ 385 and 1,155) selected from the low- 
Reo range in the experiment, for which transition was 
apparently observed. 

We observe a region in which the cylinder wake 
perturbs the boundary layer by advecting low-speed 
fluid from the wall into the outer region, and high- 
speed fluid from the core towards the wall, as rollers 
of alternating rotation are convert ed above the bound- 
ary layer. Three-dimensional disturbances appear in 
the boundary layer in correspondence to kinks in the 
spanwise rollers. The appearance of streamwise vor- 
ticity in the wake is visually well-correlated with the 
presence of low- and high-speed streaks near the wall. 
This, and the fact that these streaks are convec ted 
downstream at the freestream velocity (rather than a 
convection velocity typical of near-wall phenomena) 
indicate that the perturbations in the boundary layer 
are forced from the outside, and not generated within 
the boundary layer itself. As the wake decays, the 
advection effect becomes less significant, but the wake 


generates high levels of freestream turbulence, ffiat are 
expected to affect the flow downstream. In the low- 
Reynolds-number case the flow tends towards a more 
quiescent state towards the end of the computational 
domain and the statistics never reach turbulent levels. 
In the high-Reynolds-number case, on the other hand, 
the appearance of coherent, quasi-streamwise vortices 
is observed, which might be the prelude to the es- 
tablishment of a self- generating cycle of sweeps and 
ejections. This was supported by the establishment 
of Reynolds-stress levels comparable to those observed 
in turbulent flows. 

While some features of the computed flow fields are 
in satisfactory agreement with the measured data, the 
simulations indicate that transition does not really 
take place at the lower Reynolds number: the bound- 
ary layer is simply buffeted by the unsteady Karman 
vortex street shed by the cylinder. For Rep = 385 the 
fluctuations generated in this manner cannot be self- 
sustained due to low-Reynolds-number effects, and the 
boundary-layer flow does not transition to turbulence. 
We believe that the cause behind the discrepancy 
between the observations of Kyriakides et a l. 4 and 
the above results is to be found in an unsatisfactory 
criterion for transition detection in the experiment, 
and in the lack of sufficiently detailed measurements 
of the type particularly required for transition stud- 
ies. In contrast, in the higher Reynolds-number case, 
boundary-layer fluctuations persist even after the wake 
has decayed. This can be attributed, in part, to the 
spatial growth of the boundary layer during the in- 
teraction, which brings the Reg values closer to the 
typical range of minimum critical Reynolds numbers 
required for self-sustained turbulence in zero-pressure- 
gradient flows; some evidence could be observed that 
a self-sustaining turbulence generation cycle was be- 
ginning to be established. 

Near-term extensions to this work will involve addi- 
tional simulations to clarify the flow behavior down- 
stream of the present domain. We also plan to inves- 
tigate the possibility of higher- simulations, which 
will also allow the present study of strong wake inter- 
actions to be extended to weak interactions. It will 
be interesting to And out the extent of similarities be- 
tween the weak interactions and the familiar case of 
bypass transition due to moderate levels of fr ee-stream 
turbulence. 
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11 OF 12 


American Institute 


Aeronautics and Astro: 


'ics Paper 2003-0975 


OF 






4 Xyriakides N.K., Kastrinakis, E.G., Nychas, S.G., Goulas, A. 
1996 Proc. Inst. Mech. Eng. 210 , 167. 

5 Xyriakides N.K., Kastrinakis, E.G., Nychas, S.G., Goulas, A. 
1999 AIAA J 37 , 1197. 

6 Kyriakides N.K., Kastrinakis, E.G., Nychas, S.G., Goulas, A. 
1999 Eur J. Mech. B-Fluid Id, 1049. 

7 Chorin, A. J. 1968 Math. Comput. 22, 745. 

3 Kim, J. and Moin, ?. 1985 J. Comput . Phys. 59 308. 
9 Meneveau, C., Lund, T. S., and Cabot, W. H, 1996 J. Fluid 
Mech. 319 353-385. 

10 Mohd-Yusof, J. 1997 in CTR Armu. Res. Briefs 1997, NASA 

Ames/Stanford University, 317. 

u ?adlun, E.A., Yerrlcco, R., Criandi, ?., Mohd-Yusof, J. 2000 
J. Comput. Phys. 161 35 . 

12 Verzi cco, R., Mohd-Yusof, J., Orlandi, P., Haworth, D. 2000 
AIAA J. 38 427. 

13 Balaras, E. 2002 Submitted to Comput. and Fluids. 
14 Balaras, E., Benocci, C., and Piomelii, U. 1995 Theoret. Com- 
put. Fluid Dyn. 7 207. 

15 Balaras, E., Piomelii, U., and Wallace, J.M. 2001 J. Fluid 
Mech., 446, 1. 

16 Piomeili, U., Balaras, E., & Pascarelli, A. 2000 J. Turbulence 

1 , 001 . 

17 Orlanski, I. 1976 J. Comput Phys. 21, 251, 

18 Voke, P., and Yang, Z. 1995 Phys . Fluids 7, 2256. 

19 Roach, P.E., and Brierley, D.H. 1992 In Numerical Simulation 
of unsteady flows and transition to turbulence , O. Pironneau, W. 
Rodi, I.L. Rhyming, A.M. Savill and T.V. Truong, eds. Cam- 
bridge, 318. 

20 Wallace, J.M., Brodkey, R.S., and Sckelman, H 1972 J. Fluid 
Mech. 54, 39. 


12 OF 12 


American Institute of Aeronautics and Astronautics Paper 2003-0975 



